Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/oper_sys/weyl/weyl_lsp.lha/lisp-numbers.lisp

Text File | 1991-10-29 | 12.6 KB | 493 lines

;;; -*- Mode:Lisp; Package:Weyli; Base:10; Lowercase:T; Syntax:Common-Lisp -*-
;;; ===========================================================================
;;;                  Lisp Numbers
;;; ===========================================================================
;;; (c) Copyright 1989, 1991 Cornell University

;;; $Id: lisp-numbers.lisp,v 2.21 1991/10/30 00:00:39 rz Exp $


(in-package "WEYLI")

(define-domain-element-classes floating-point-numbers float)

(defmethod print-object ((d floating-point-numbers) stream)
  #-Genera
  (format stream "R")
  #+Genera
  (format stream "~'bR~"))

(defvar *floating-point-numbers* (make-instance 'floating-point-numbers))

;; This isn't quite right yet, since I don't have the special type of
;; floating point number yet.
;; The class: REAL-NUMBERS is defined in bigfloat.lisp
(define-domain-creator real-numbers (&optional (precision 16))
  (if (= precision 16)
      *floating-point-numbers*
      (make-instance 'real-numbers :precision precision)))

;; This is slightly inefficient, but who cares...  I want to localize
;; the knowledge of how to create domains in the MAKE-...* functions.
(defun get-real-numbers (&optional (precision 16))
  (if (= precision 16)
      (add-domain (lambda (d)
            (eql d *floating-point-numbers*))
    (make-real-numbers* precision))
      (add-domain (lambda (d)
        (and (eql (class-name (class-of d)) 'real-numbers)
             (eql (fp-precision d) precision)))
    (make-real-numbers* precision))))

;; The Lisp number domain is a special field whose elements are
;; represented as Lisp numbers.  This class is unique.

(define-domain-element-classes lisp-numbers
    integer rational)

(defmethod print-object ((d lisp-numbers) stream)
  #+Genera
  (format stream "~'bQlisp~")
  #-Genera
  (princ "Qlisp"  stream))

#+IGNORE  ; Defined in domain-support...
(defvar *lisp-numbers* (make-instance 'lisp-numbers))

(defun get-lisp-numbers ()
  (or *lisp-numbers*
      (error "You need to RESET-DOMAINS!")))

(defun make-lisp-numbers ()
  (get-lisp-numbers))

(defmethod domain-of ((element float))
  *floating-point-numbers*)

(defmethod domain-of ((element (or integer ratio lisp:complex)))
  *lisp-numbers*)

(defmethod number? ((element t))
  nil)

(defmethod number? ((element number))
  t)

;; Special printer for complex numbers to that they are more readable.
(defmethod print-object ((c lisp:complex) stream)
  (cond ((lisp:zerop (imagpart c))
     (princ (realpart c) stream ))
    ((lisp:zerop (realpart c))
     (format stream "(~S i)" (imagpart c)))
    (t (format stream "(~S + ~S i)" (realpart c) (imagpart c)))))

;;For documentation purposes  (provided by Lisp)
#+ignore
(defmethod print-object ((n number) stream)
  (princ n stream))

(defmethod numerator ((r integer))
  r)

(defmethod denominator ((r integer))
  1)

(defmethod numerator ((r ratio))
  (lisp::numerator r))

(defmethod denominator ((r ratio))
  (lisp::denominator r))

(defmethod < ((x number) (y number))
  (lisp:< x y))

(defmethod <= ((x number) (y number))
  (lisp:<= x y))

(defmethod = ((x number) (y number))
  (lisp:= x y))

(defmethod >= ((x number) (y number))
  (lisp:>= x y))

(defmethod > ((x number) (y number))
  (lisp:> x y))

(defmethod max-pair ((x number) (y number))
  (if (lisp:> x y) x y))

(defmethod min-pair ((x number) (y number))
  (if (lisp:> x y) y x))

(defmethod 0? ((x number))
  (lisp:zerop x))

(defmethod 1? ((x number))
  (= x 1))

(defmethod zero ((r lisp-numbers))
  0)

(defmethod one ((r lisp-numbers))
  1)

(defmethod minus ((x number))
  (lisp:- x))

(defmethod minus? ((x number))
  (lisp:minusp x))

(defmethod minus? ((x lisp:complex))
  nil)

(defmethod plus? ((x number))
  (lisp:plusp x))

(defmethod plus? ((x lisp:complex))
  (not (lisp:zerop x)))

(defmethod abs ((x number))
  (lisp:abs x))

;; The explicit mention of integer is needed to overide the coercions
;; in coercions.lisp

(defmethod plus ((x (or integer number)) (y (or integer number)))
  (lisp:+ x y))

(defmethod difference ((x (or integer number)) (y (or integer number)))
  (lisp:- x y))

(defmethod times ((x (or integer number)) (y (or integer number)))
  (lisp:* x y))

(defmethod recip ((x number))
  (lisp:/ x))

(defmethod expt ((n number) (e number))
  (lisp:expt n e))

(defmethod quotient ((a (or integer number)) (b (or integer number)))
  (lisp:/ a b))

(defmethod floor ((a number) &optional b)
  (if b
      (lisp:floor a b)
      (lisp:floor a)))

(defmethod ceiling ((a number) &optional b)
  (if b
      (lisp:ceiling a b)
      (lisp:ceiling a)))

(defmethod truncate ((a number) &optional b)
  (if b
      (lisp:truncate a b)
      (lisp:truncate a)))

(defmethod round ((a number) &optional b)
  (if b
      (lisp:round a b)
      (lisp:round  a)))

(defmethod remainder ((a number) (b number))
  (lisp:rem a b))

(defmethod gcd ((a integer) (b integer))
  (lisp:gcd a b))

;; Do we really need this???
(defmethod gcd ((a float) (b float))
  a)

(defmethod lcm ((a integer) (b integer))
  (lisp:* (lisp:/ a (lisp:gcd a b)) b))

(defun faster-isqrt (n)
  (let (n-len-quarter n-half n-half-isqrt init-value q r iterated-value)
    "argument n must be a non-negative integer"
    (cond
      ((> n 24)
       ;; theoretically (> n 7) ,i.e., n-len-quarter > 0
       (setq n-len-quarter (ash (integer-length n) -2))
       (setq n-half (ash n (- (ash n-len-quarter 1))))
       (setq n-half-isqrt (faster-isqrt n-half))
       (setq init-value (ash (1+ n-half-isqrt) n-len-quarter))
       (multiple-value-setq (q r) (floor n init-value))
       (setq iterated-value (ash (+ init-value q) -1))
       (if (eq (logbitp 0 q) (logbitp 0 init-value)) ; same sign
       ;; average is exact and we need to test the result
       (let ((m (- iterated-value init-value)))
         (if (> (* m m) r)
         (- iterated-value 1)
         iterated-value))
       ;; average was not exact, we take value
       iterated-value))
      ((> n 15) 4)
      ((> n  8) 3)
      ((> n  3) 2)
      ((> n  0) 1)
      ((> n -1) 0)
      (t nil))))

(defvar *big-primes* ()
  "List of large primes by decending size")

;; Return the next prime less than its argument, and that fits into a
;; word.
(defun newprime (p)
  (do ((pl *big-primes* (cdr pl)))
      ((null pl) (setq p (find-smaller-prime p))
       (setq *big-primes* (nconc *big-primes* (list p)))
       p)
    (if (lisp:< (car pl) p) (return (car pl)))))

;; Rabin's probabilistic primality algorithm isn't used here because it
;; isn't much faster than the simple one for numbers about the size of a
;; word.
(defun find-smaller-prime (p)
  "Finds biggest prime less than fixnum p"
  (if (evenp p) (setq p (1- p)))
  (loop for pp = (lisp:- p 2) then (lisp:- pp 2) until (lisp:< pp 0)
    when (prime? pp)
      do (return pp)))

(defun repeated-squaring (mult one)
  (lambda (base exp)
    (if (lisp:zerop exp) one
    (let ((prod one))
      (loop
        (if (oddp exp)
        (setq prod (%funcall mult prod base)))
        (setq exp (truncate exp 2))
        (if (lisp:zerop exp)
        (return prod))
        (setq base (%funcall mult base base)))))))

(defmethod power-of? ((m integer) &optional n)
  (cond ((typep n 'integer)
     (loop for test = n then (* test n)
           for i upfrom 1
           when (= test m)
         do (return (values n i))
           when (> test m)
         do (return nil)))
    (t (error "Haven't implemented the rest of the cases"))))

;; These two really should be in GFP, but because of LUCID braindamage,
;; they have to be here to avoid warnings.

(defun reduce-modulo-integer (value modulus)
  (unless (lisp:zerop modulus)
    (setq value (lisp:rem value modulus)))
  (if (lisp:< value 0) (lisp:+ value modulus)
      value))

(defun expt-modulo-integer (base expt modulus)
  (%funcall (repeated-squaring
          (lambda (a b) (reduce-modulo-integer (lisp:* a b) modulus))
          1)
        base expt))

(defmethod prime? ((p integer))
  (and (or (lisp:< p 14.)
       (and (lisp:= 1 (expt-modulo-integer 13. (1- p) p))
        (lisp:= 1 (expt-modulo-integer 3 (1- p) p))))
       (null (cdr (setq p (factor p))))
       (lisp:= 1 (cdar p))))

(defun all-divisors (n)
  (let ((factors (factor n)))
    (loop with divisors = (list 1)
      for (prime . times) in factors
      do (loop for i from 1 to times
           appending (loop for divisor in divisors
                   collect (* divisor (lisp:expt prime i)))
             into temp
           finally (setq divisors (append temp divisors)))
         finally (return (sort divisors #'lisp:<)))))

(defvar *factor-method* 'simple-integer-factor)

(defmacro count-multiple-integer-factors (N divisor)
  `(loop with i = 0
     do (multiple-value-bind (quo rem) (lisp:truncate ,N ,divisor)
          (when (not (lisp:zerop rem))
        (if (not (lisp:zerop i))
            (push (cons ,divisor i) ans))
        (return t))
          (setq ,N quo)
          (incf i))))

(defmethod factor ((N integer))
  (let ((*factor-method* *factor-method*)
    ans factors)
    (when (lisp:minusp N)
      (push (cons -1 1) ans)
      (setq N (lisp:- N)))
    (count-multiple-integer-factors N 2)
    (count-multiple-integer-factors N 3)
    (count-multiple-integer-factors N 5)
    (unless (lisp::= N 1)
      (loop
    (multiple-value-setq (N factors) (%funcall *factor-method* N))
    (setq ans (append factors ans))
    (if (lisp::= N 1) (return t))))
    (uniformize-factor-list ans)))

(defun uniformize-factor-list (ans)
  (loop for pairs on (sort ans (lambda (a b) (< (first a) (first b))))
    when (or (null (rest pairs))
         (not (lisp::= (first (first pairs))
             (first (second pairs)))))
      collect (first pairs)
    else do (incf (rest (second pairs)))))

;; In general each factorization method should return just one factor.

(defvar *skip-chain-for-3-and-5* (circular-list 4 2 4 2 4 6 2 6))

(defun simple-integer-factor (N)
  (let ((increments *skip-chain-for-3-and-5*)
    (divisor 7)
    ans)
    (flet ((simple-integer-factor-internal (N)
         (let ((limit (lisp:isqrt N)))
           (loop
         (cond ((lisp:= N 1)
            (return (values N ans)))
               ((lisp:> divisor limit)
            (return (values 1 (list (cons N 1)))))
               (t (count-multiple-integer-factors N divisor)))
         (setq divisor (lisp:+ divisor (pop increments)))))))
      (setq *factor-method* #'simple-integer-factor-internal)
      (simple-integer-factor-internal N))))

(defun fermat-integer-factor (N)
  (loop for x = (1+ (lisp:isqrt N)) then (+ x 1)
    for w = (lisp:- (lisp:* x x) N)
    for y = (lisp:isqrt w)
    do (when (lisp:zerop (lisp:- w (lisp:* y y)))
         (let ((u (lisp:+ x y))
           (v (lisp:- x y)))
           (return (if (1? v)
               (values 1 (list (cons u 1)))
               (values u (factor v))))))))

#| Knuth's addition-subtraction version of Fermat's algorithm |

(defun fermat-integer-factor (N)
  (let* ((isqrt-N (lisp:isqrt N))
     (x (1+ (* 2 isqrt-N)))
     (y 1)
     (r (- (* isqrt-N isqrt-N) N)))
    (loop
      (cond ((= r 0)
         (return
           (let ((f (/ (+ x y -2) 2))
             (g (/ (- x y) 2)))
         (if (= g 1)
             (values 1 (list (cons f 1)))
             (values 1 (append (factor f) (factor g)))))))
        ((< r 0)
         (incf r x)
         (incf x 2)))
      (decf r y)
      (incf y 2))))

(defun list-of-primes (N)
  (cons 2
    (loop for p upfrom 3 by 2 below N
          when (prime? p) collect p)))

(defun make-integer-GCD-list (max-prime size-limit)
  (let ((GCD-list ()))
    (loop for p in (list-of-primes max-prime)
      with prod = 1 and prime-list = ()
      do (setq prod (* prod p))
         (cond ((> prod size-limit)
            (push (list (/ prod p) prime-list)
              GCD-list)
            (setq prod p)
            (setq prime-list (list p)))
           (t (push p prime-list))))
    GCD-list))


||#

(defun totient (x)
  (do ((factors (factor x) (rest factors))
       (totient 1 (lisp:* totient
              (lisp:- (lisp:expt (caar factors) (cdar factors))
                  (lisp:expt (caar factors) (1- (cdar factors)))))))
      ((null factors)
       totient)))

(defmethod sqrt ((x number))
  (lisp:sqrt  x))

(defmethod sin ((x number))
  (lisp:sin x))

(defmethod cos ((x number))
  (lisp:cos x))

(defmethod tan ((x number))
  (lisp:tan x))

(defmethod asin ((x number))
  (lisp:asin x))

(defmethod acos ((x number))
  (lisp:acos x))

(defmethod atan ((x number) &optional y)
  (cond ((null y)
     (lisp:atan x))
    ((numberp y)
     (lisp:atan x y))
    (t (atan (coerce x (domain-of y)) y))))

(defmethod sinh ((x number))
  (lisp:sinh x))

(defmethod cosh ((x number))
  (lisp:cosh x))

(defmethod tanh ((x number))
  (tanh x))

(defmethod asinh ((x number))
  (lisp:asinh x))

(defmethod acosh ((x number))
  (lisp:acosh x))

(defmethod atanh ((x number))
  (lisp:atanh x))

(defmethod exp ((x number))
  (lisp:exp x))

(defmethod log2 ((x number) (base number))
  (lisp:log x base))

(defmethod log ((x number))
  (lisp:log x))

(defmethod conjugate  ((x number))
  (lisp:conjugate x))

(defmethod realpart ((x number))
  (lisp:realpart x))

(defmethod imagpart ((x number))
  (lisp:imagpart x))

(defmethod phase ((x number))
  (lisp:phase x))

(defmethod signum ((x number))
  (lisp:signum x))